library(tidyverse) ; library(reshape2) ; library(glue) ; library(plotly) ; library(dendextend)
library(RColorBrewer) ; library(viridis) ; require(gridExtra) ; library(GGally) ; library(ggpubr)
library(expss)
library(polycor)
library(foreach) ; library(doParallel)
library(knitr) ; library(kableExtra)
library(biomaRt)
library(clusterProfiler) ; library(ReactomePA) ; library(DOSE) ; library(org.Hs.eg.db)
library(WGCNA)

SFARI_colour_hue = function(r) {
  pal = c('#FF7631','#FFB100','#E8E328','#8CC83F','#62CCA6','#59B9C9','#b3b3b3','#808080','gray','#d9d9d9')[r]
}

Load preprocessed dataset (preprocessing code in 01_data_preprocessing.Rmd) and clustering (pipeline in 05_WGCNA.Rmd)

# Gandal dataset
load('./../Data/preprocessed_data.RData')
datExpr = datExpr %>% data.frame
rownames(datExpr) = datGenes$ensembl_gene_id
datMeta = datMeta %>% mutate(ID = paste0('X',description))


# GO Neuronal annotations: regex 'neuron' in GO functional annotations and label the genes that make a match as neuronal
GO_annotations = read.csv('./../Data/genes_GO_annotations.csv')
GO_neuronal = GO_annotations %>% filter(grepl('neuron', go_term)) %>% 
              mutate('ID'=as.character(ensembl_gene_id)) %>% 
              dplyr::select(-ensembl_gene_id) %>% distinct(ID) %>%
              mutate('Neuronal'=1)


# SFARI Genes
SFARI_genes = read_csv('./../../../SFARI/Data/SFARI_genes_01-03-2020_w_ensembl_IDs.csv')
SFARI_genes = SFARI_genes[!duplicated(SFARI_genes$ID) & !is.na(SFARI_genes$ID),]

# Old SFARI Genes
old_SFARI_genes = read_csv('./../../../SFARI/Data/SFARI_genes_08-29-2019_w_ensembl_IDs.csv')
old_SFARI_genes = old_SFARI_genes[!duplicated(old_SFARI_genes$ID) & !is.na(old_SFARI_genes$ID),]


# Clusterings
clusterings = read_csv('./../Data/clusters.csv')


# Update DE_info with SFARI and Neuronal information
genes_info = DE_info %>% mutate('ID'=datGenes$ensembl_gene_id, padj = adj.P.Val, log2FoldChange = logFC) %>% 
             left_join(SFARI_genes, by='ID') %>% 
             mutate(`gene-score`=ifelse(is.na(`gene-score`), 'Others', `gene-score`)) %>%
             left_join(old_SFARI_genes %>% rename(`gene-score` = 'old-gene-score') %>% 
                       dplyr::select(ID, `old-gene-score`), by = 'ID') %>%
             mutate(`old-gene-score`=ifelse(is.na(`old-gene-score`), 'Others', `old-gene-score`)) %>%
             left_join(GO_neuronal, by='ID') %>% left_join(clusterings, by='ID') %>%
             mutate(Neuronal=ifelse(is.na(Neuronal), 0, Neuronal)) %>%
             mutate(gene.score=ifelse(`gene-score`=='Others' & Neuronal==1, 'Neuronal', `gene-score`), 
                    old.gene.score=ifelse(`old-gene-score`=='Others' & Neuronal==1,'Neuronal',`old-gene-score`), 
                    significant=padj<0.05 & !is.na(padj))

# Add gene symbol
getinfo = c('ensembl_gene_id','external_gene_id')
mart = useMart(biomart='ENSEMBL_MART_ENSEMBL', dataset='hsapiens_gene_ensembl',
               host='feb2014.archive.ensembl.org') ## Gencode v19
gene_names = getBM(attributes=getinfo, filters=c('ensembl_gene_id'), values=genes_info$ID, mart=mart)

genes_info = genes_info %>% left_join(gene_names, by=c('ID'='ensembl_gene_id'))


clustering_selected = 'DynamicHybrid'
genes_info$Module = genes_info[,clustering_selected]

dataset = read.csv(paste0('./../Data/dataset_', clustering_selected, '.csv'))
dataset$Module = dataset[,clustering_selected]
dataset$gene.score = as.character(dataset$gene.score)
dataset$gene.score[dataset$gene.score=='None'] = 'Others' 


rm(DE_info, GO_annotations, clusterings, getinfo, mart, efit)


Selecting Top Modules


plot_data = dataset %>% dplyr::select(Module, MTcor) %>% distinct %>% 
            mutate(alpha = ifelse(abs(MTcor)>0.85, 1, 0.3))

top_modules = plot_data %>% arrange(desc(MTcor)) %>% filter(abs(MTcor)>0.85) %>% pull(Module) %>% as.character

Selecting Modules with a Module-Diagnosis correlation magnitude larger than 0.85 (instead of 0.9 as we did with Gandal’s dataset because the relation is not as strong in this dataset)

The 2 modules that fulfill this condition are #79AF00, #00B813

ggplotly(plot_data %>% ggplot(aes(reorder(Module, -MTcor), MTcor)) + 
         geom_bar(stat='identity', fill = plot_data$Module, alpha = plot_data$alpha) + 
         geom_hline(yintercept =c(0.85, -0.85), color = 'gray', linetype = 'dotted') + 
         xlab('Modules')+ ylab('Module-Diagnosis Correlation') + theme_minimal() + 
         theme(axis.text.x = element_text(angle = 90, hjust = 1)))


The modules consist mainly of points with high (absolute) values in PC2 (which we know is related to LFC), so this result is consistent with the high correlation between Module and Diagnosis, although some of the points with the highest PC2 values do not belong to these top modules

The genes belonging to the modules with the positive Module-Diagnosis correlation have positive LFC values and the genes belonging to the modules with the negative Module-Diagnosis correlation have negative values

pca = datExpr %>% prcomp

plot_data = data.frame('ID'=rownames(datExpr), 'PC1' = pca$x[,1], 'PC2' = pca$x[,2]) %>%
            left_join(dataset, by='ID') %>% 
            left_join(genes_info %>% dplyr::select(ID, external_gene_id), by='ID') %>%
            dplyr::select(ID, external_gene_id, PC1, PC2, Module, gene.score) %>%
            mutate(ImportantModules = ifelse(Module %in% top_modules, as.character(Module), 'Others')) %>%
            mutate(color = ifelse(ImportantModules=='Others','gray',ImportantModules),
                   alpha = ifelse(ImportantModules=='Others', 0.1, 0.6),
                   gene_id = paste0(ID, ' (', external_gene_id, ')')) %>%
            apply_labels(ImportantModules = 'Top Modules')

cro(plot_data$ImportantModules)
 #Total 
 Top Modules 
   #00B813  708
   #79AF00  399
   Others  14351
   #Total cases  15458
ggplotly(plot_data %>% ggplot(aes(PC1, PC2, color=ImportantModules)) + 
         geom_point(alpha=plot_data$alpha, color=plot_data$color, aes(ID=gene_id)) + theme_minimal() +
         xlab(paste0('PC1 (',round(100*summary(pca)$importance[2,1],2),'%)')) +
         ylab(paste0('PC2 (',round(100*summary(pca)$importance[2,2],2),'%)')) +
         ggtitle('Genes belonging to the Modules with the strongest relation to ASD'))
rm(pca)




SFARI Genes


List of top 20 SFARI Genes in top modules ordered by SFARI score and Gene Significance

list_top_SFARI_genes = function(table_data, module) {
  
  t = table_data %>% filter(Module == module & `SFARI score` %in% c(1,2,3)) %>% 
      slice_head(n=20) %>% dplyr::select(-Module, -`Ensembl ID`)
 
  return(t)
}


table_data = dataset %>% left_join(genes_info %>% dplyr::select(ID, external_gene_id), by='ID') %>%
             dplyr::select(ID, external_gene_id, GS, gene.score, Module) %>% 
             arrange(gene.score, desc(abs(GS))) %>%
             dplyr::rename('Ensembl ID' = ID, 'Gene Symbol' = external_gene_id, 
                    'SFARI score' = gene.score, 'Gene Significance' = GS)

kable(list_top_SFARI_genes(table_data, top_modules[1]), 
      caption=paste0('Top SFARI Genes for Module ', top_modules[1])) %>%
      kable_styling(full_width = F)
Top SFARI Genes for Module #79AF00
Gene Symbol Gene Significance SFARI score
CASK 0.7991471 1
DLG4 0.6492071 1
SRCAP 0.3073226 1
MECP2 0.0496312 1
DEAF1 0.2414273 2
TAF6 0.2175760 2
GABRG3 0.1697446 2
OXTR 0.1531516 2
CDH13 0.1529404 2
SNX5 0.0422536 2
CMIP 0.4673635 3
MUC4 0.4548544 3
MBD3 0.4369314 3
AVPR1B 0.4227676 3
GRM5 0.3342928 3
PSD3 0.3328927 3
HTR3A 0.2924446 3
LRRC1 0.2344388 3
ZC3H11A 0.2164467 3
TTN 0.1221977 3
kable(list_top_SFARI_genes(table_data, top_modules[2]), 
      caption=paste0('Top SFARI Genes for Module ', top_modules[2])) %>%
      kable_styling(full_width = F)
Top SFARI Genes for Module #00B813
Gene Symbol Gene Significance SFARI score
SON -0.3279955 1
SHANK2 -0.2358767 1
ARX -0.2314617 1
VPS13B -0.2273199 1
PTPN11 -0.1303289 1
KDM3B -0.1218145 1
ICA1 -0.6947087 2
GGNBP2 -0.6446193 2
CACNA2D3 -0.6297103 2
TERF2 -0.5561551 2
MTOR -0.4602444 2
RBFOX1 -0.4470894 2
KCNS3 -0.4198392 2
GALNT8 -0.2812489 2
APBB1 -0.2747579 2
CHRNA7 -0.2721684 2
VIL1 -0.2406558 2
SLC7A3 -0.2062360 2
OTUD7A -0.1204609 2
APBA2 -0.6365906 3
rm(table_data, list_top_SFARI_genes)




Module Eigengenes


Since these modules have the strongest relation to autism, this pattern should be reflected in their model eigengenes, having two different behaviours for the samples corresponding to autism and the ones corresponding to control

In all cases, the Eigengenes separate the behaviour between autism and control samples very clearly (p-value < \(10^{-4}\)). Modules with positive Module-Diagnosis correlation have higher eigengenes in the ASD samples and Modules with a negative correlation, in the Control samples

plot_EGs = function(module){

  plot_data = data.frame('ID' = rownames(MEs), 'MEs' = MEs[,paste0('ME', module)], 
                         'Diagnosis' = datMeta$Diagnosis)

  p = plot_data %>% ggplot(aes(Diagnosis, MEs, fill=Diagnosis)) + 
      geom_boxplot(outlier.colour='#cccccc', outlier.shape='o', outlier.size=3) +
      ggtitle(paste0('Module ', module, '  (MTcor=',round(dataset$MTcor[dataset$Module == module][1],2),')')) + 
      stat_compare_means(method = 't.test', method.args = list(var.equal = FALSE), label = 'p.signif',
                        ref.group = 'ASD') +
      ylab('Module Eigengenes') + theme_minimal() + theme(legend.position='none')
  
  return(p)
}


# Calculate Module Eigengenes
ME_object = datExpr %>% t %>% moduleEigengenes(colors = genes_info$Module)
MEs = orderMEs(ME_object$eigengenes)

p1 = plot_EGs(top_modules[1])
p2 = plot_EGs(top_modules[2])

grid.arrange(p1, p2, nrow=1)

rm(plot_EGs, ME_object, MEs, p1, p2)




Identifying representative genes for each Module


In the WGCNA pipeline, the most representative genes of each module are selected based on having a high module membership as well as a high (absolute) gene significance, so I’m going to do the same

SFARI Genes don’t seem to be more representative than the rest of the genes

create_plot = function(module){
  
  plot_data = dataset %>% dplyr::select(ID, paste0('MM.',gsub('#','',module)), GS, gene.score) %>% 
              filter(dataset$Module==module)
  colnames(plot_data)[2] = 'Module'
  
  SFARI_colors = as.numeric(names(table(as.character(plot_data$gene.score)[plot_data$gene.score!='Others'])))
  
  p = ggplotly(plot_data %>% mutate(gene.score = ifelse(gene.score =='Others', 'Not in SFARI', gene.score)) %>% 
               ggplot(aes(Module, GS, color=gene.score)) +
               geom_point(alpha=0.5, aes(ID=ID)) +  xlab('Module Membership') + ylab('Gene Significance') +
               ggtitle(paste0('Module ', module, '  (MTcor = ', 
                              round(dataset$MTcor[dataset$Module == module][1],2),')')) +
               scale_color_manual(values=SFARI_colour_hue(r=c(SFARI_colors,7))) + theme_minimal())
  
  return(p)
}

create_plot(top_modules[1])
create_plot(top_modules[2])
rm(create_plot)

Top 20 genes for each module


Ordered by \(\frac{MM+|GS|}{2}\)

There aren’t that many SFARI genes in the top genes of the modules

create_table = function(module){
  top_genes = dataset %>% left_join(genes_info %>% dplyr::select(ID, external_gene_id), by='ID') %>% 
              dplyr::select(ID, external_gene_id, paste0('MM.',gsub('#','',module)), GS, gene.score) %>% 
              filter(dataset$Module==module) %>% dplyr::rename('MM' = paste0('MM.',gsub('#','',module))) %>% 
              mutate(Relevance = (MM+abs(GS))/2, 
                     gene.score = ifelse(gene.score =='Others', 'Not in SFARI', gene.score)) %>% 
              arrange(by=-Relevance) %>% top_n(20) %>% 
              dplyr::rename('Gene Symbol' = external_gene_id, 'SFARI Score' = gene.score)
  return(top_genes)
}

top_genes = list()
for(i in 1:length(top_modules)) top_genes[[i]] = create_table(top_modules[i])

kable(top_genes[[1]] %>% dplyr::select(-ID), caption=paste0('Top 20 genes for Module ', top_modules[1], 
      '  (MTcor = ', round(dataset$MTcor[dataset$Module == top_modules[1]][1],2),')')) %>% 
      kable_styling(full_width = F)
Top 20 genes for Module #79AF00 (MTcor = 0.87)
Gene Symbol MM GS SFARI Score Relevance
WWC1 0.7770464 0.8108212 Not in SFARI 0.7939338
CASK 0.7886111 0.7991471 1 0.7938791
CCDC50 0.8122779 0.7306045 Not in SFARI 0.7714412
ABLIM3 0.7420932 0.7196454 Not in SFARI 0.7308693
RIMS4 0.7437251 0.6842736 Not in SFARI 0.7139994
DLG4 0.7218259 0.6492071 1 0.6855165
RAB11FIP4 0.6854802 0.6449110 Not in SFARI 0.6651956
JUN 0.6378685 0.6890684 Not in SFARI 0.6634684
PRR5 0.6026740 0.7175114 Not in SFARI 0.6600927
MXD4 0.6532816 0.6281044 Not in SFARI 0.6406930
RBM39 0.7523953 0.5196978 Not in SFARI 0.6360465
CCDC102A 0.6342699 0.6362624 Not in SFARI 0.6352661
SERTAD1 0.6300838 0.6313406 Not in SFARI 0.6307122
TMEM132A 0.6166367 0.6436432 Not in SFARI 0.6301399
TNFRSF1A 0.7071841 0.5462540 Not in SFARI 0.6267190
SBF2 0.6532819 0.5634664 Not in SFARI 0.6083741
KIAA1522 0.5645191 0.6477544 Not in SFARI 0.6061367
COMMD2 0.5308295 0.6683542 Not in SFARI 0.5995919
KCNK12 0.6686929 0.5259704 Not in SFARI 0.5973317
MSN 0.6179045 0.5757439 Not in SFARI 0.5968242
kable(top_genes[[2]] %>% dplyr::select(-ID), caption=paste0('Top 20 genes for Module ', top_modules[2], 
      '  (MTcor = ', round(dataset$MTcor[dataset$Module == top_modules[2]][1],2),')')) %>% 
      kable_styling(full_width = F)
Top 20 genes for Module #00B813 (MTcor = -0.91)
Gene Symbol MM GS SFARI Score Relevance
CIRBP 0.8238270 -0.8158820 Not in SFARI 0.8198545
EIF3K 0.7973822 -0.7788495 Not in SFARI 0.7881159
FIBP 0.7705627 -0.8052324 Not in SFARI 0.7878976
UQCRC1 0.8000285 -0.7699530 Not in SFARI 0.7849907
COPG1 0.7859792 -0.7813419 Not in SFARI 0.7836606
PFKP 0.7868055 -0.7437745 Not in SFARI 0.7652900
RTN4IP1 0.8020323 -0.7283473 Not in SFARI 0.7651898
EIF2AK1 0.7984542 -0.7119941 Not in SFARI 0.7552241
PYGB 0.7675212 -0.7128820 Not in SFARI 0.7402016
C17orf49 0.7426747 -0.7237342 Not in SFARI 0.7332044
ICA1 0.7663316 -0.6947087 2 0.7305202
TCEA2 0.7690770 -0.6881259 Not in SFARI 0.7286015
MRPL2 0.7539501 -0.7031087 Not in SFARI 0.7285294
ACTR1B 0.6882590 -0.7683051 Not in SFARI 0.7282821
ALOX12B 0.6670812 -0.7739453 Not in SFARI 0.7205133
OGDHL 0.7048090 -0.7346210 Not in SFARI 0.7197150
PRC1 0.7844209 -0.6470912 Not in SFARI 0.7157560
GLS2 0.7180064 -0.7038931 Not in SFARI 0.7109497
FH 0.6930093 -0.7237622 Not in SFARI 0.7083858
IMMT 0.6590283 -0.7498505 Not in SFARI 0.7044394
rm(create_table, i)
pca = datExpr %>% prcomp

ids = c()
for(tg in top_genes) ids = c(ids, tg$ID)

plot_data = data.frame('ID'=rownames(datExpr), 'PC1' = pca$x[,1], 'PC2' = pca$x[,2]) %>%
            left_join(dataset, by='ID') %>% dplyr::select(ID, PC1, PC2, Module, gene.score) %>%
            mutate(color = ifelse(Module %in% top_modules, as.character(Module), 'gray')) %>%
            mutate(alpha = ifelse(color %in% top_modules & ID %in% ids, 1, 0.1))

plot_data %>% ggplot(aes(PC1, PC2)) + geom_point(alpha=plot_data$alpha, color=plot_data$color) + 
              xlab(paste0('PC1 (',round(100*summary(pca)$importance[2,1],2),'%)')) +
              ylab(paste0('PC2 (',round(100*summary(pca)$importance[2,2],2),'%)')) +
              theme_minimal() + ggtitle('Most relevant genes for top Modules')

rm(ids, pca, tg, plot_data)

Level of expression by Diagnosis for top genes, ordered by relevance (defined above): There is a visible difference in level of expression between diagnosis groups in all of these genes

create_plot = function(i){
  
  plot_data = datExpr[rownames(datExpr) %in% top_genes[[i]]$ID,] %>% mutate('ID' = rownames(.)) %>% 
              melt(id.vars='ID') %>% mutate(variable = gsub('X','',variable)) %>%
              left_join(top_genes[[i]], by='ID') %>%
              left_join(datMeta %>% dplyr::select(description, Diagnosis),
                        by = c('variable'='description')) %>% arrange(desc(Relevance))
  
  p = ggplotly(plot_data %>% mutate(external_gene_id=factor(`Gene Symbol`, 
                                    levels=unique(plot_data$`Gene Symbol`), ordered=T)) %>%
               ggplot(aes(`Gene Symbol`, value, fill=Diagnosis)) + geom_boxplot() + theme_minimal() +
                      ggtitle(paste0('Top Genes for Module ', top_modules[i], ' (MTcor = ',
                      round(dataset$MTcor[dataset$Module == top_modules[i]][1],2), ')')) + 
                      xlab('') + ylab('Level of Expression') +
                      theme(axis.text.x = element_text(angle = 90, hjust = 1)))
  return(p)
}

create_plot(1)
create_plot(2)
rm(create_plot)




Enrichment Analysis


Using the package clusterProfiler. Performing Gene Set Enrichment Analysis (GSEA) and Over Representation Analysis (ORA) using the following datasets:

file_name = './../Data/GSEA.RData'

if(file.exists(file_name)){
  load(file_name)
} else {
  ##############################################################################
  # PREPARE DATASET
  
  # Create dataset with top modules membership and removing the genes without an assigned module
  EA_dataset = data.frame('ensembl_gene_id' = genes_info$ID, module = genes_info$Module)  %>% 
               filter(genes_info$Module!='gray')
  
  # Assign Entrez Gene Id to each gene
  getinfo = c('ensembl_gene_id','entrezgene')
  mart = useMart(biomart='ENSEMBL_MART_ENSEMBL', dataset='hsapiens_gene_ensembl', 
                 host='feb2014.archive.ensembl.org')
  biomart_output = getBM(attributes=getinfo, filters=c('ensembl_gene_id'), 
                         values=genes_info$ID[genes_info$Module!='gray'], mart=mart)
  
  EA_dataset = biomart_output %>% dplyr::rename('ID' = ensembl_gene_id) %>%
               left_join(dataset %>% dplyr::select(ID, contains('MM.')), by='ID')

  
  ##############################################################################
  # PERFORM ENRICHMENT
  
  # Following https://yulab-smu.github.io/clusterProfiler-book/chapter8.html
  
  modules = dataset$Module[dataset$Module!='gray'] %>% as.character %>% table %>% names
  nPerm = 1e5 # 100 times more than the default
  
  enrichment_GO = list()         # Gene Ontology
  enrichment_DO = list()         # Disease Ontology
  enrichment_DGN = list()        # Disease Gene Networks
  enrichment_KEGG = list()       # Kyoto Encyclopedia of Genes and Genomes
  enrichment_Reactome = list()   # Reactome: Pathway db
  
  
  for(module in modules){
    cat('\n')
    cat(paste0('Module: ', which(modules == module), '/', length(modules)))
    geneList = EA_dataset[,paste0('MM.',substring(module,2))]
    names(geneList) = EA_dataset[,'entrezgene'] %>% as.character
    geneList = sort(geneList, decreasing = TRUE)
    
    enrichment_GO[[module]] = gseGO(geneList, OrgDb = org.Hs.eg.db, pAdjustMethod = 'bonferroni', 
                                    pvalueCutoff = 0.1, nPerm = nPerm, verbose = FALSE, seed = TRUE) %>% 
                              data.frame
    enrichment_DO[[module]] = gseDO(geneList, pAdjustMethod = 'bonferroni', pvalueCutoff = 0.1,
                                    nPerm = nPerm, verbose = FALSE, seed = TRUE) %>% data.frame
    enrichment_DGN[[module]] = gseDGN(geneList, pAdjustMethod = 'bonferroni', pvalueCutoff = 0.1,
                                      nPerm = nPerm, verbose = FALSE, seed = TRUE) %>% data.frame
    enrichment_KEGG[[module]] = gseKEGG(geneList, organism = 'human', pAdjustMethod = 'bonferroni', 
                                        pvalueCutoff = 0.1, nPerm = nPerm, verbose = FALSE, seed = TRUE) %>% 
                                data.frame
    enrichment_Reactome[[module]] = gsePathway(geneList, organism = 'human', pAdjustMethod = 'bonferroni', 
                                               pvalueCutoff = 0.1, nPerm = nPerm, verbose = F, seed = T) %>% 
                                    data.frame
    
    # Temporal save, just in case SFARI Genes enrichment fails
    save(enrichment_GO, enrichment_DO, enrichment_DGN, enrichment_KEGG, enrichment_Reactome, file=file_name)
  }
  

  ##############################################################################
  # PERFROM ENRICHMENT FOR SFARI GENES
  
  # BUILD MAPPING BETWEEN GENES AND SFARI

  # Build TERM2GENE: dataframe of 2 columns with term and gene
  term2gene = biomart_output %>% 
              left_join(genes_info %>% dplyr::select(ID, `gene-score`), 
                         by = c('ensembl_gene_id'='ID')) %>% dplyr::select(-ensembl_gene_id) %>% 
              mutate('SFARI' = ifelse(`gene-score` != 'Others','SFARI','Others'),
                     `gene-score` = ifelse(`gene-score` != 'Others', 
                                           paste0('SFARI Score ',`gene-score`), 'Others')) %>%
              melt(id.vars = 'entrezgene') %>% dplyr::select(value, entrezgene) %>% 
              dplyr::rename('term' = value, 'gene' = entrezgene) %>% distinct
  
  
  # PERFORM GSEA
  enrichment_SFARI = list()
  cat('\n\nGSEA OF SFARI GENES\n')
  
  for(module in modules){
    cat('\n')
    cat(paste0('Module: ', which(modules == module), '/', length(modules)))
    geneList = EA_dataset[,paste0('MM.',substring(module,2))]
    names(geneList) = EA_dataset[,'entrezgene'] %>% as.character
    geneList = sort(geneList, decreasing = TRUE)
      
    enrichment_SFARI[[module]] = clusterProfiler::GSEA(geneList, pAdjustMethod = 'bonferroni',  nPerm = nPerm,
                                                       TERM2GENE = term2gene, pvalueCutoff=1, maxGSSize=2e3,
                                                        verbose = FALSE, seed = TRUE) %>% data.frame
    
    # Temporal save
    save(enrichment_GO, enrichment_DO, enrichment_DGN, enrichment_KEGG, enrichment_Reactome, 
         enrichment_SFARI, file=file_name)
  }
  
  ##############################################################################
  # PERFROM ENRICHMENT FOR OLD SFARI GENES
  
  # BUILD MAPPING BETWEEN GENES AND SFARI

  # Build TERM2GENE: dataframe of 2 columns with term and gene
  term2gene = biomart_output %>% 
              left_join(genes_info %>% dplyr::select(ID, `old-gene-score`), 
                         by = c('ensembl_gene_id'='ID')) %>% dplyr::select(-ensembl_gene_id) %>% 
              mutate('SFARI' = ifelse(`old-gene-score` != 'Others','SFARI','Others'),
                     `old-gene-score` = ifelse(`old-gene-score` != 'Others', 
                                           paste0('SFARI Score ',`old-gene-score`), 'Others')) %>%
              melt(id.vars = 'entrezgene') %>% dplyr::select(value, entrezgene) %>% 
              dplyr::rename('term' = value, 'gene' = entrezgene) %>% distinct
  
  
  # PERFORM GSEA
  enrichment_old_SFARI = list()
  cat('\n\nGSEA OF OLD SFARI GENES\n')
  
  for(module in modules){
    cat('\n')
    cat(paste0('Module: ', which(modules == module), '/', length(modules)))
    geneList = EA_dataset[,paste0('MM.',substring(module,2))]
    names(geneList) = EA_dataset[,'entrezgene'] %>% as.character
    geneList = sort(geneList, decreasing = TRUE)
      
    enrichment_old_SFARI[[module]] = clusterProfiler::GSEA(geneList, pAdjustMethod = 'bonferroni',  
                                                           nPerm = nPerm, TERM2GENE = term2gene, pvalueCutoff=1, 
                                                           maxGSSize=2e3, verbose = FALSE, seed = TRUE) %>% 
                                     data.frame
    
    # Temporal save
    save(enrichment_GO, enrichment_DO, enrichment_DGN, enrichment_KEGG, enrichment_Reactome, 
         enrichment_SFARI, enrichment_old_SFARI, file=file_name)
  }
  ##############################################################################
  # Save enrichment results
  save(enrichment_GO, enrichment_DO, enrichment_DGN, enrichment_KEGG, enrichment_Reactome, 
       enrichment_SFARI, enrichment_old_SFARI, file=file_name)
  
  rm(getinfo, mart, biomart_output, module, term2gene, geneList, EA_dataset, nPerm)
}

# Rename lists
GSEA_GO = enrichment_GO
GSEA_DGN = enrichment_DGN
GSEA_DO = enrichment_DO
GSEA_KEGG = enrichment_KEGG
GSEA_Reactome = enrichment_Reactome
GSEA_SFARI = enrichment_SFARI
GSEA_old_SFARI = enrichment_old_SFARI
file_name = './../Data/ORA.RData'

if(file.exists(file_name)){
  load(file_name)
} else {
  
  ##############################################################################
  # PREPARE DATASET
  
  # Create dataset with top modules membership and removing the genes without an assigned module
  EA_dataset = data.frame('ensembl_gene_id' = genes_info$ID, module = genes_info$Module)  %>% 
               filter(genes_info$Module!='gray')
  
  # Assign Entrez Gene Id to each gene
  getinfo = c('ensembl_gene_id','entrezgene')
  mart = useMart(biomart='ENSEMBL_MART_ENSEMBL', dataset='hsapiens_gene_ensembl', 
                 host='feb2014.archive.ensembl.org')
  biomart_output = getBM(attributes=getinfo, filters=c('ensembl_gene_id'), 
                         values=genes_info$ID[genes_info$Module!='gray'], mart=mart)
  
  EA_dataset = biomart_output %>% dplyr::rename('ID' = ensembl_gene_id) %>%
               left_join(dataset %>% dplyr::select(ID, Module), by='ID')

  
  ##############################################################################
  # PERFORM ENRICHMENT
  
  # Following https://yulab-smu.github.io/clusterProfiler-book/chapter8.html
  
  modules = dataset$Module[dataset$Module!='gray'] %>% as.character %>% table %>% names
  universe = EA_dataset$entrezgene %>% as.character
  
  enrichment_GO = list()         # Gene Ontology
  enrichment_DO = list()         # Disease Ontology
  enrichment_DGN = list()        # Disease Gene Networks
  enrichment_KEGG = list()       # Kyoto Encyclopedia of Genes and Genomes
  enrichment_Reactome = list()   # Reactome: Pathway db
  
  
  for(module in modules){
    
    genes_in_module = EA_dataset$entrezgene[EA_dataset$Module == module]
    
    enrichment_GO[[module]] = enrichGO(gene = genes_in_module, universe = universe, OrgDb = org.Hs.eg.db, 
                                       ont = 'All', pAdjustMethod = 'bonferroni', pvalueCutoff = 0.1, 
                                       qvalueCutoff = 1) %>% data.frame
    enrichment_DO[[module]] = enrichDO(gene = genes_in_module, universe = universe, qvalueCutoff = 1,
                                       pAdjustMethod = 'bonferroni', pvalueCutoff = 0.1) %>% data.frame
    enrichment_DGN[[module]] = enrichDGN(gene = genes_in_module, universe = universe, qvalueCutoff = 1,
                                         pAdjustMethod = 'bonferroni', pvalueCutoff = 0.1) %>% data.frame
    enrichment_KEGG[[module]] = enrichKEGG(gene = genes_in_module, universe = universe, qvalueCutoff = 1,
                                           pAdjustMethod = 'bonferroni', pvalueCutoff = 0.1) %>% data.frame
    enrichment_Reactome[[module]] = enrichPathway(gene = genes_in_module, universe = universe, qvalueCutoff = 1,
                                                  pAdjustMethod = 'bonferroni', pvalueCutoff = 0.1) %>% 
                                    data.frame
  }
  
  # Temporal save, just in case SFARI Genes enrichment fails
  save(enrichment_GO, enrichment_DO, enrichment_DGN, enrichment_KEGG, enrichment_Reactome, file=file_name)
  
  
  ##############################################################################
  # PERFROM ENRICHMENT FOR SFARI GENES
  
  # BUILD MAPPING BETWEEN GENES AND SFARI

  # Build TERM2GENE: dataframe of 2 columns with term and gene
  term2gene = biomart_output %>% 
              left_join(genes_info %>% dplyr::select(ID, `gene-score`), 
                         by = c('ensembl_gene_id'='ID')) %>% dplyr::select(-ensembl_gene_id) %>% 
              mutate('SFARI' = ifelse(`gene-score` != 'Others','SFARI','Others'),
                     `gene-score` = ifelse(`gene-score` != 'Others', 
                                           paste0('SFARI Score ',`gene-score`), 'Others')) %>%
              melt(id.vars = 'entrezgene') %>% dplyr::select(value, entrezgene) %>% 
              dplyr::rename('term' = value, 'gene' = entrezgene) %>% distinct
  
  
  # PERFORM GSEA
  enrichment_SFARI = list()
  
  for(module in modules){
      genes_in_module = EA_dataset$entrezgene[EA_dataset$Module == module]
      
      enrichment_SFARI[[module]] = enricher(gene = genes_in_module, universe = universe, 
                                            pAdjustMethod = 'bonferroni', TERM2GENE = term2gene, 
                                            pvalueCutoff = 1, qvalueCutoff = 1, maxGSSize = 50000) %>% 
                                    data.frame %>% dplyr::select(-geneID,-Description)
  }

  ##############################################################################
  # PERFROM ENRICHMENT FOR SFARI GENES
  
  # BUILD MAPPING BETWEEN GENES AND SFARI

  # Build TERM2GENE: dataframe of 2 columns with term and gene
  term2gene = biomart_output %>% 
              left_join(genes_info %>% dplyr::select(ID, `old-gene-score`), 
                         by = c('ensembl_gene_id'='ID')) %>% dplyr::select(-ensembl_gene_id) %>% 
              mutate('SFARI' = ifelse(`old-gene-score` != 'Others','SFARI','Others'),
                     `old-gene-score` = ifelse(`old-gene-score` != 'Others', 
                                           paste0('SFARI Score ',`old-gene-score`), 'Others')) %>%
              melt(id.vars = 'entrezgene') %>% dplyr::select(value, entrezgene) %>% 
              dplyr::rename('term' = value, 'gene' = entrezgene) %>% distinct
  
  
  # PERFORM GSEA
  enrichment_old_SFARI = list()
  
  for(module in modules){
      genes_in_module = EA_dataset$entrezgene[EA_dataset$Module == module]
      
      enrichment_old_SFARI[[module]] = enricher(gene = genes_in_module, universe = universe, 
                                                pAdjustMethod = 'bonferroni', TERM2GENE = term2gene, 
                                                pvalueCutoff = 1, qvalueCutoff = 1, maxGSSize = 5e4) %>% 
                                       data.frame %>% dplyr::select(-geneID,-Description)
  }
  
  ##############################################################################
  # Save enrichment results
  save(enrichment_GO, enrichment_DO, enrichment_DGN, enrichment_KEGG, enrichment_Reactome, 
       enrichment_SFARI, enrichment_old_SFARI, file=file_name)
  
  rm(getinfo, mart, biomart_output, gene, module, term2gene, genes_in_module, EA_dataset, universe, file_name)
}

# Rename lists
ORA_GO = enrichment_GO
ORA_DGN = enrichment_DGN
ORA_DO = enrichment_DO
ORA_KEGG = enrichment_KEGG
ORA_Reactome = enrichment_Reactome
ORA_SFARI = enrichment_SFARI
ORA_old_SFARI = enrichment_old_SFARI

rm(enrichment_GO, enrichment_DO, enrichment_DGN, enrichment_KEGG, enrichment_Reactome, enrichment_SFARI, 
   enrichment_old_SFARI)


Both GSEA and ORA are commonly used to study enrichment in sets of genes, but when using them for studying our modules both have shortcomings:

modules = dataset$Module[dataset$Module!='gray'] %>% as.character %>% table %>% names

module = sample(modules,1)

plot_data = dataset %>% dplyr::select(Module, paste0('MM.',gsub('#','',module))) %>% 
            mutate(in_module = substring(Module,2) == gsub('#','',module), selected_module = module) %>%
            mutate(alpha = ifelse(in_module, 0.8, 0.1))
colnames(plot_data)[2] = 'MM'

p = plot_data %>% ggplot(aes(selected_module, MM, color = in_module)) + geom_jitter(alpha = plot_data$alpha) + 
    xlab('') + ylab('Module Membership') + coord_flip() + theme_minimal() + 
    theme(legend.position = 'none')

ggExtra::ggMarginal(p, type = 'density', groupColour = TRUE, groupFill = TRUE, margins = 'x', size=1)

rm(modules, module, p, plot_data)

So perhaps it could be useful to use both methods together, since they seem to complement each other’s shortcomings very well, performing the enrichment using both methods and identifying the terms that are found to be enriched by both

Note: Since the enrichment in both methods is quite a stric restriction, we decide to relax the corrected p-value threshold (using Bonferroni correction) to 0.1.

compare_methods = function(GSEA_list, ORA_list){
  
  for(top_module in top_modules){
  
    cat(paste0('  \n  \nEnrichments for Module ', top_module, ' (MTcor=',
               round(dataset$MTcor[dataset$Module==top_module][1],2), '):  \n  \n'))
    
    GSEA = GSEA_list[[top_module]]
    ORA = ORA_list[[top_module]]
    
    cat(paste0('GSEA has ', nrow(GSEA), ' enriched terms  \n'))
    cat(paste0('ORA has  ', nrow(ORA), ' enriched terms  \n'))
    cat(paste0(sum(ORA$ID %in% GSEA$ID), ' terms are enriched in both methods  \n  \n'))

    plot_data = GSEA %>% mutate(pval_GSEA = p.adjust) %>% dplyr::select(ID, Description, NES, pval_GSEA) %>%
                inner_join(ORA %>% mutate(pval_ORA = p.adjust) %>% 
                           dplyr::select(ID, pval_ORA, GeneRatio, qvalue), by = 'ID') 
    
    if(nrow(plot_data)>0){
      print(plot_data %>% mutate(pval_mean = pval_ORA + pval_GSEA) %>% 
                          arrange(pval_mean) %>% dplyr::select(-pval_mean) %>% 
            kable %>% kable_styling(full_width = F))
    }
  } 
}


plot_results = function(GSEA_list, ORA_list){
  
  l = htmltools::tagList()

  for(i in 1:length(top_modules)){
    
    GSEA = GSEA_list[[top_modules[i]]]
    ORA = ORA_list[[top_modules[i]]]
    
    plot_data = GSEA %>% mutate(pval_GSEA = p.adjust) %>% dplyr::select(ID, Description, NES, pval_GSEA) %>%
                inner_join(ORA %>% mutate(pval_ORA = p.adjust) %>% dplyr::select(ID, pval_ORA), by = 'ID')
    
    if(nrow(plot_data)>5){
      min_val = min(min(plot_data$pval_GSEA), min(plot_data$pval_ORA))
      max_val = max(max(max(plot_data$pval_GSEA), max(plot_data$pval_ORA)),0.05)
      ggp = ggplotly(plot_data %>% ggplot(aes(pval_GSEA, pval_ORA, color = NES)) + 
                     geom_point(aes(id = Description)) + 
                     geom_vline(xintercept = 0.05, color = 'gray', linetype = 'dotted') + 
                     geom_hline(yintercept = 0.05, color = 'gray', linetype = 'dotted') + 
                     ggtitle(paste0('Enriched terms in common for Module ', top_modules[i])) +
                     scale_x_continuous(limits = c(min_val, max_val)) + 
                     scale_y_continuous(limits = c(min_val, max_val)) + 
                     xlab('Corrected p-value for GSEA') + ylab('Corrected p-value for ORA') +
                     scale_colour_viridis(direction = -1) + theme_minimal() + coord_fixed())
      l[[i]] = ggp
    }
  }
  
  return(l)
}


KEGG

compare_methods(GSEA_KEGG, ORA_KEGG)

Enrichments for Module #79AF00 (MTcor=0.87):

GSEA has 34 enriched terms
ORA has 3 enriched terms
3 terms are enriched in both methods

ID Description NES pval_GSEA pval_ORA GeneRatio qvalue
hsa05322 Systemic lupus erythematosus 2.383671 0.0076646 0.0000729 15/196 3.48e-05
hsa05203 Viral carcinogenesis 2.039923 0.0081770 0.0000695 20/196 3.48e-05
hsa05034 Alcoholism 2.512285 0.0080476 0.0002834 18/196 9.02e-05

Enrichments for Module #00B813 (MTcor=-0.91):

GSEA has 18 enriched terms
ORA has 0 enriched terms
0 terms are enriched in both methods


Reactome

compare_methods(GSEA_Reactome, ORA_Reactome)

Enrichments for Module #79AF00 (MTcor=0.87):

GSEA has 94 enriched terms
ORA has 73 enriched terms
64 terms are enriched in both methods

ID Description NES pval_GSEA pval_ORA GeneRatio qvalue
R-HSA-171306 Packaging Of Telomere Ends 2.797288 0.0288531 0.0000000 14/248 0.0000000
R-HSA-110330 Recognition and association of DNA glycosylase with site containing an affected purine 2.809859 0.0289917 0.0000000 14/248 0.0000000
R-HSA-110331 Cleavage of the damaged purine 2.809859 0.0289917 0.0000000 14/248 0.0000000
R-HSA-73927 Depurination 2.809859 0.0289917 0.0000000 14/248 0.0000000
R-HSA-110328 Recognition and association of DNA glycosylase with site containing an affected pyrimidine 2.690112 0.0291747 0.0000000 14/248 0.0000000
R-HSA-110329 Cleavage of the damaged pyrimidine 2.690112 0.0291747 0.0000000 14/248 0.0000000
R-HSA-73928 Depyrimidination 2.690112 0.0291747 0.0000000 14/248 0.0000000
R-HSA-3214842 HDMs demethylate histones 2.459638 0.0292524 0.0000000 16/248 0.0000000
R-HSA-73929 Base-Excision Repair, AP Site Formation 2.684661 0.0292524 0.0000000 14/248 0.0000000
R-HSA-73728 RNA Polymerase I Promoter Opening 2.946047 0.0293359 0.0000000 14/248 0.0000000
R-HSA-5334118 DNA methylation 2.918911 0.0294813 0.0000001 14/248 0.0000000
R-HSA-427359 SIRT1 negatively regulates rRNA expression 2.759346 0.0295267 0.0000001 14/248 0.0000000
R-HSA-5625886 Activated PKN1 stimulates transcription of AR (androgen receptor) regulated genes KLK2 and KLK3 2.863767 0.0295749 0.0000002 14/248 0.0000000
R-HSA-5693571 Nonhomologous End-Joining (NHEJ) 2.303682 0.0295919 0.0000000 15/248 0.0000000
R-HSA-606279 Deposition of new CENPA-containing nucleosomes at the centromere 2.682341 0.0295919 0.0000002 14/248 0.0000000
R-HSA-774815 Nucleosome assembly 2.682341 0.0295919 0.0000002 14/248 0.0000000
R-HSA-1221632 Meiotic synapsis 2.611073 0.0296705 0.0000003 14/248 0.0000000
R-HSA-2299718 Condensation of Prophase Chromosomes 2.945113 0.0297640 0.0000005 14/248 0.0000000
R-HSA-212300 PRC2 methylates histones and DNA 2.784749 0.0297640 0.0000005 14/248 0.0000000
R-HSA-2559586 DNA Damage/Telomere Stress Induced Senescence 2.500765 0.0298261 0.0000001 15/248 0.0000000
R-HSA-5693565 Recruitment and ATM-mediated phosphorylation of repair and signaling proteins at DNA double strand breaks 2.401013 0.0298393 0.0000001 15/248 0.0000000
R-HSA-157579 Telomere Maintenance 2.304934 0.0298393 0.0000014 14/248 0.0000000
R-HSA-5693606 DNA Double Strand Break Response 2.400206 0.0298997 0.0000002 15/248 0.0000000
R-HSA-427389 ERCC6 (CSB) and EHMT2 (G9a) positively regulate rRNA expression 2.815317 0.0298997 0.0000002 15/248 0.0000000
R-HSA-912446 Meiotic recombination 2.545407 0.0300142 0.0000001 16/248 0.0000000
R-HSA-3214858 RMTs methylate histone arginines 2.283929 0.0304344 0.0000000 17/248 0.0000000
R-HSA-73772 RNA Polymerase I Promoter Escape 2.484752 0.0304344 0.0000185 14/248 0.0000004
R-HSA-5250924 B-WICH complex positively regulates rRNA expression 2.530772 0.0304344 0.0000185 14/248 0.0000004
R-HSA-1912408 Pre-NOTCH Transcription and Translation 2.398050 0.0304712 0.0000003 16/248 0.0000000
R-HSA-201722 Formation of the beta-catenin:TCF transactivating complex 2.398704 0.0304712 0.0000026 15/248 0.0000001
R-HSA-8936459 RUNX1 regulates genes involved in megakaryocyte differentiation and platelet function 2.299876 0.0305792 0.0000039 15/248 0.0000001
R-HSA-5625740 RHO GTPases activate PKNs 2.429326 0.0306756 0.0000069 15/248 0.0000002
R-HSA-3214815 HDACs deacetylate histones 2.660197 0.0306822 0.0000006 16/248 0.0000000
R-HSA-977225 Amyloid fiber formation 2.458757 0.0308546 0.0000021 16/248 0.0000001
R-HSA-1912422 Pre-NOTCH Expression and Processing 2.327395 0.0309407 0.0000052 16/248 0.0000001
R-HSA-1500620 Meiosis 2.194344 0.0309407 0.0000052 16/248 0.0000001
R-HSA-5578749 Transcriptional regulation by small RNAs 2.122131 0.0308383 0.0001460 14/248 0.0000025
R-HSA-5250913 Positive epigenetic regulation of rRNA expression 2.450169 0.0309586 0.0000331 15/248 0.0000006
R-HSA-73854 RNA Polymerase I Promoter Clearance 2.387812 0.0309990 0.0000456 15/248 0.0000009
R-HSA-427413 NoRC negatively regulates rRNA expression 2.362381 0.0308939 0.0001704 14/248 0.0000028
R-HSA-2559582 Senescence-Associated Secretory Phenotype (SASP) 2.574914 0.0311734 0.0000011 17/248 0.0000000
R-HSA-73864 RNA Polymerase I Transcription 2.386248 0.0311425 0.0000534 15/248 0.0000010
R-HSA-5250941 Negative epigenetic regulation of rRNA expression 2.314879 0.0309407 0.0002675 14/248 0.0000044
R-HSA-5617472 Activation of anterior HOX genes in hindbrain development during early embryogenesis 2.300994 0.0312880 0.0000167 16/248 0.0000004
R-HSA-5619507 Activation of HOX genes during differentiation 2.300994 0.0312880 0.0000167 16/248 0.0000004
R-HSA-2559580 Oxidative Stress Induced Senescence 2.402468 0.0314089 0.0000308 16/248 0.0000006
R-HSA-1474165 Reproduction 2.096910 0.0314604 0.0000552 16/248 0.0000010
R-HSA-211000 Gene Silencing by RNA 1.888874 0.0315485 0.0003883 15/248 0.0000060
R-HSA-9018519 Estrogen-dependent gene expression 2.322324 0.0319318 0.0000882 17/248 0.0000016
R-HSA-3214847 HATs acetylate histones 2.025256 0.0319366 0.0003048 16/248 0.0000049
R-HSA-68875 Mitotic Prophase 2.185369 0.0319366 0.0016379 15/248 0.0000234
R-HSA-8878171 Transcriptional regulation by RUNX1 1.808700 0.0340641 0.0003342 21/248 0.0000053
R-HSA-212165 Epigenetic regulation of gene expression 2.129608 0.0319318 0.0025336 15/248 0.0000350
R-HSA-2559583 Cellular Senescence 2.065143 0.0331844 0.0014131 18/248 0.0000205
R-HSA-8939211 ESR-mediated signaling 2.398064 0.0337253 0.0017596 19/248 0.0000247
R-HSA-9006931 Signaling by Nuclear Receptors 2.100650 0.0346960 0.0251258 19/248 0.0003138
R-HSA-73884 Base Excision Repair 2.101587 0.0603541 0.0000103 14/248 0.0000002
R-HSA-69473 G2/M DNA damage checkpoint 2.059731 0.0609424 0.0000026 15/248 0.0000001
R-HSA-73886 Chromosome Maintenance 2.028023 0.0616677 0.0001063 14/248 0.0000018
R-HSA-8939236 RUNX1 regulates transcription of genes involved in differentiation of HSCs 1.859668 0.0629809 0.0004414 15/248 0.0000066
R-HSA-157118 Signaling by NOTCH 1.763084 0.0683002 0.0001010 22/248 0.0000018
R-HSA-3247509 Chromatin modifying enzymes 1.678713 0.0696326 0.0133438 20/248 0.0001739
R-HSA-4839726 Chromatin organization 1.678713 0.0696326 0.0133438 20/248 0.0001739
R-HSA-3214841 PKMTs methylate histone lysines 1.986687 0.0897702 0.0000029 14/248 0.0000001

Enrichments for Module #00B813 (MTcor=-0.91):

GSEA has 27 enriched terms
ORA has 0 enriched terms
0 terms are enriched in both methods


Gene Ontology

compare_methods(GSEA_GO, ORA_GO)

Enrichments for Module #79AF00 (MTcor=0.87):

GSEA has 31 enriched terms
ORA has 65 enriched terms
0 terms are enriched in both methods

Enrichments for Module #00B813 (MTcor=-0.91):

GSEA has 29 enriched terms
ORA has 1 enriched terms
0 terms are enriched in both methods


Disease Ontology

compare_methods(GSEA_DO, ORA_DO)

Enrichments for Module #79AF00 (MTcor=0.87):

GSEA has 29 enriched terms
ORA has 0 enriched terms
0 terms are enriched in both methods

Enrichments for Module #00B813 (MTcor=-0.91):

GSEA has 18 enriched terms
ORA has 0 enriched terms
0 terms are enriched in both methods


Disease Gene Network

compare_methods(GSEA_DGN, ORA_DGN)

Enrichments for Module #79AF00 (MTcor=0.87):

GSEA has 6 enriched terms
ORA has 8 enriched terms
4 terms are enriched in both methods

ID Description NES pval_GSEA pval_ORA GeneRatio qvalue
umls:C0278504 Non-small cell lung cancer stage I 2.110432 0.0732876 0.00e+00 16/346 0.0e+00
umls:C0278484 Malignant neoplasm of colon stage IV 2.180269 0.0739020 1.00e-07 16/346 0.0e+00
umls:C0206659 Embryonal Carcinoma 1.955250 0.0789265 3.00e-06 20/346 6.0e-07
umls:C0036421 Systemic Scleroderma 1.668825 0.0973237 5.65e-05 36/346 8.8e-06

Enrichments for Module #00B813 (MTcor=-0.91):

GSEA has 0 enriched terms
ORA has 0 enriched terms
0 terms are enriched in both methods



Session info

sessionInfo()
## R version 3.6.3 (2020-02-29)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.4 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
## 
## locale:
##  [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_GB.UTF-8        LC_COLLATE=en_GB.UTF-8    
##  [5] LC_MONETARY=en_GB.UTF-8    LC_MESSAGES=en_GB.UTF-8   
##  [7] LC_PAPER=en_GB.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats4    parallel  stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] WGCNA_1.69             fastcluster_1.1.25     dynamicTreeCut_1.63-1 
##  [4] org.Hs.eg.db_3.8.2     AnnotationDbi_1.46.1   IRanges_2.18.3        
##  [7] S4Vectors_0.22.1       Biobase_2.44.0         BiocGenerics_0.30.0   
## [10] DOSE_3.10.2            ReactomePA_1.28.0      clusterProfiler_3.12.0
## [13] biomaRt_2.40.5         kableExtra_1.1.0       knitr_1.28            
## [16] doParallel_1.0.15      iterators_1.0.12       foreach_1.5.0         
## [19] polycor_0.7-10         expss_0.10.2           ggpubr_0.2.5          
## [22] magrittr_1.5           GGally_1.5.0           gridExtra_2.3         
## [25] viridis_0.5.1          viridisLite_0.3.0      RColorBrewer_1.1-2    
## [28] dendextend_1.13.4      plotly_4.9.2           glue_1.4.1            
## [31] reshape2_1.4.4         forcats_0.5.0          stringr_1.4.0         
## [34] dplyr_1.0.0            purrr_0.3.4            readr_1.3.1           
## [37] tidyr_1.1.0            tibble_3.0.1           ggplot2_3.3.2         
## [40] tidyverse_1.3.0       
## 
## loaded via a namespace (and not attached):
##   [1] tidyselect_1.1.0      RSQLite_2.2.0         htmlwidgets_1.5.1    
##   [4] grid_3.6.3            BiocParallel_1.18.1   munsell_0.5.0        
##   [7] codetools_0.2-16      preprocessCore_1.46.0 miniUI_0.1.1.1       
##  [10] withr_2.2.0           colorspace_1.4-1      GOSemSim_2.10.0      
##  [13] highr_0.8             rstudioapi_0.11       ggsignif_0.6.0       
##  [16] labeling_0.3          urltools_1.7.3        polyclip_1.10-0      
##  [19] bit64_0.9-7           farver_2.0.3          vctrs_0.3.1          
##  [22] generics_0.0.2        xfun_0.12             R6_2.4.1             
##  [25] graphlayouts_0.7.0    bitops_1.0-6          reshape_0.8.8        
##  [28] fgsea_1.10.1          gridGraphics_0.5-0    assertthat_0.2.1     
##  [31] promises_1.1.0        scales_1.1.1          ggraph_2.0.3         
##  [34] nnet_7.3-14           enrichplot_1.4.0      ggExtra_0.9          
##  [37] gtable_0.3.0          tidygraph_1.2.0       rlang_0.4.6          
##  [40] splines_3.6.3         lazyeval_0.2.2        acepack_1.4.1        
##  [43] impute_1.58.0         broom_0.5.5           europepmc_0.4        
##  [46] checkmate_2.0.0       BiocManager_1.30.10   yaml_2.2.1           
##  [49] modelr_0.1.6          crosstalk_1.1.0.1     backports_1.1.8      
##  [52] httpuv_1.5.2          qvalue_2.16.0         Hmisc_4.4-0          
##  [55] tools_3.6.3           ggplotify_0.0.5       ellipsis_0.3.1       
##  [58] ggridges_0.5.2        Rcpp_1.0.4.6          plyr_1.8.6           
##  [61] base64enc_0.1-3       progress_1.2.2        RCurl_1.98-1.2       
##  [64] prettyunits_1.1.1     rpart_4.1-15          cowplot_1.0.0        
##  [67] haven_2.2.0           ggrepel_0.8.2         cluster_2.1.0        
##  [70] fs_1.4.0              data.table_1.12.8     DO.db_2.9            
##  [73] triebeard_0.3.0       reprex_0.3.0          reactome.db_1.68.0   
##  [76] matrixStats_0.56.0    xtable_1.8-4          mime_0.9             
##  [79] hms_0.5.3             evaluate_0.14         XML_3.99-0.3         
##  [82] jpeg_0.1-8.1          readxl_1.3.1          compiler_3.6.3       
##  [85] crayon_1.3.4          htmltools_0.4.0       later_1.0.0          
##  [88] Formula_1.2-3         lubridate_1.7.4       DBI_1.1.0            
##  [91] tweenr_1.0.1          dbplyr_1.4.2          MASS_7.3-51.6        
##  [94] rappdirs_0.3.1        Matrix_1.2-18         cli_2.0.2            
##  [97] igraph_1.2.5          pkgconfig_2.0.3       rvcheck_0.1.8        
## [100] foreign_0.8-76        xml2_1.2.5            webshot_0.5.2        
## [103] rvest_0.3.5           digest_0.6.25         graph_1.62.0         
## [106] rmarkdown_2.1         cellranger_1.1.0      fastmatch_1.1-0      
## [109] htmlTable_1.13.3      curl_4.3              shiny_1.4.0.2        
## [112] graphite_1.30.0       lifecycle_0.2.0       nlme_3.1-147         
## [115] jsonlite_1.7.0        fansi_0.4.1           pillar_1.4.4         
## [118] lattice_0.20-41       fastmap_1.0.1         httr_1.4.1           
## [121] survival_3.1-12       GO.db_3.8.2           UpSetR_1.4.0         
## [124] png_0.1-7             bit_1.1-15.2          ggforce_0.3.1        
## [127] stringi_1.4.6         blob_1.2.1            latticeExtra_0.6-29  
## [130] memoise_1.1.0